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Beyond chemical accuracy: The pseudopotential approximation in diffusion Monte 
Carlo calculations of the HCP to BCC phase transition in beryllium 
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Motivated by the disagreement between recent diffusion Monte Carlo calculations and experiments 
on the phase transition pressure between the ambient and beta-Sn phases of silicon, we present a 
study of the HCP to BCC phase transition in beryllium. This lighter element provides an oppor¬ 
tunity for directly testing many of the approximations required for calculations on silicon and may 
suggest a path towards increasing the practical accuracy of diffusion Monte Carlo calculations of 
solids in general. We demonstrate that the single largest approximation in these calculations is the 
pseudopotential approximation. After removing this we find excellent agreement with experiment 
for the ambient HCP phase and results similar to careful calculations using density functional theory 
for the phase transition pressure. 
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While diffusion quantum Monte Carlo (DMC) has 
proven to be among the most accurate methods for elec¬ 
tronic structure calculations^P— some recent results have 
shown that the method may have some difficulty in pre¬ 
dicting energetics between competing structures in con¬ 
densed phase materials heavier than hydrogen.— Al¬ 
though the method is systematically improvable, these 
studies have all employed two important approximations. 
The first of these is the use of pseudopotentials to remove 
the need to explicitly treat the chemically inert core elec¬ 
trons, and the second is use of a single Slater-determinant 
trial wavefunction to restrict the nodal surface that con¬ 
strains the diffusion Monte Carlo calculations. Given the 
computational expense of DMC, it is crucial for the fu¬ 
ture use of the method to understand the relative impor¬ 
tance of these approximations and if possible understand 
how to reduce their magnitude. 

In this paper, we address the importance of these ap¬ 
proximations for the conceptually simple but computa¬ 
tionally challenging solid-solid phase transition of beryl¬ 
lium from the ambient HCP phase to BCC at high pres¬ 
sure. Beryllium is an interesting material due to its im¬ 
portance as a lightweight material with relevance to a 
number of high pressure applications, including the re¬ 
cently proposed MagLIF and NIF fusion concepts— De¬ 
spite this technological significance, the phase diagram 
of beryllium at high temperatures and pressures is not 
tightly constrained from experiments despite a great deal 
of effort— The material has also attracted numerous 
attempts to calculate its phase diagram using electronic 
structure techniques.— In addition, the phase transi¬ 
tion from HCP to BCC structures is extraordinarily sen¬ 
sitive to any errors that may be present in the underlying 
electronic structure calculations. 

To illustrate this, we present a series of density func¬ 
tional theory (DFT) based calculations of the zero tem¬ 
perature phase transition pressure. The DFT framework 
requires an approximation in the choice of a functional 
to treat the effects of electron correlation. Although it 
is possible to classify the complexity of approximations 
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FIG. 1. (color online) Calculated zero temperature phase 
transition pressure from HCP to BCC. Calculations were per¬ 
formed using the LDA functional and include a quasiharmonic 
treatment of the zero point motion of the beryllium. The dot 
labeled LDA denotes the calculated energy difference between 
the zero pressure volumes of BCC and HCP. The line shows 
the change in transition pressure as this difference is artifi- 
cally changed. The values of the energy difference between the 
equilibrium density HCP and BCC phases and corresponding 
transition pressures are also plotted for DMC with different 
choices for treating the core electrons. 

to this exchange correlation functional)^ increasing this 
complexity does not guarantee increased accuracy. In 
fact, there is no a priori way to select the most accurate 
functional for a specific calculation and the spread in pre¬ 
dictions of the total energy in a solid can often be larger 
than the 1 kcal/mol (~ 43 meV or 0.8 mHa / atom) 
target that is typically deemed chemical accuracy. Here 
both the HCP and BCC phases are treated within the 
local density approximation, utilizing the quasiharmonic 
approximation to account for the zero-point motion of 
the relatively light nuclei. This initial calculation results 
in a phase transition pressure of ss 410 GPa, shown figure 
[T] as the point labeled LDA. However, we also consider 
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the simplest possible error between the phases: a rigid 
shift of the energies of one phase compared to another, 
shown by the line in the same figure. This produces an 
enormous change in the phase transition pressure. In 
fact, as this energy difference is shifted by 0.3 kcal/mol 
(less than 1/3 of the value typically desired for chemical 
accuracy) the phase transition pressure changes by over 
150 GPa. 

This extreme sensitivity motivates application of dif¬ 
fusion Monte Carlo (DMC) to the problem. While the 
size of the errors in DMC calculations of real materials 
has yet to be seen, the approximation made in DMC is of 
an entirely different class (being primarily a topological 
constraint on the many body wavefunction) than the one 
made in density functional theory. For this reason one 
may seek by DMC to build a diversity of approximations, 
easing the uncertainty in predictions caused by reliance 
on a single computational technique. 

As a starting point for these DMC calculations, we 
utilize the methodology presented in our recent paper 
on applying quantum Monte Carlo to solidsii Summa¬ 
rizing this briefly, we use the QMCPACK code to perform 
DMC calculations on periodic systems.— Pseudopo¬ 
tentials are generated using the opium pseudopotential 
generation progrant^^ so as to accurately reproduce the 
results of all electron LAPW calculations for the lat¬ 
tice constant and bulk modulus of the ambient phase of 
the material. For DMC, the pseudopotential was evalu¬ 
ated using a variational technique— One and two-body 
Jastrow factors were used to construct a Slater-Jastrow 
wavefunction where the single particle orbitals were ex¬ 
tracted from DFT calculations using the local density 
approximation. The timestep and the number of super¬ 
cell twists were chosen so that the total energy of the 
calculations was converged to within 1 meV per atom. 
To accurately calculate a phase transition pressure, the 
two body finite size effects need to be converged abso¬ 
lutely rather than relatively within each phase. In the 
case of HCP beryllium, this required 132 atom supercells 
with 2-body finite size errors estimated using a combina¬ 
tion of the model periodic Coulomb interaction^^ and a 
kinetic energy correction estimated from the long range 
piece of the variational Monte Carlo optimized 2-body 
Jastrow factor— Any remaining error in the calculations 
will be due either to the fixed node approximation or to 
the construction and evaluation of the pseudopotential. 

The c/a ratio of HCP beryllium must be optimized for 
each volume. This was done optimizing the DMC energy 
over volume conserving strains, yielding c/a ratios shown 
in Fig 121 The calculated c/a ratio at ambient pressure is 
1.569 -b/- 0.004, which is in excellent agreement with the 
experimentally measured 1.568— 

We now calculate the energy as a function of volume of 
HCP beryllium for the optimized structures, fitting the 
sum of resulting values and the quasiharmonic zero point 
energy determined within the LDA to a Vinet equation 
of state— We find an equilibrium volume that is slightly 
too small at 52.27 ± 0.02 bohr^/atom compared to the 
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FIG. 2. (color online) Calculated HCP c/a from quadratic 
fits to the DMC energy vs volume conserving strain is shown 
in red with error bars. Experimental c/a versus specific vol¬ 
ume taken from diamond anvil cell experiments reported in 
reference [TJ is plotted as blue circles. 


experimental 54.776 bohr^/atom. The bulk modulus is 
somewhat 1 too large at 124.21 ± 0.37 CPa, compared to 
the experimental value of 116.8 CPa. Despite this disap¬ 
pointing performance, we continue calculating the same 
energy versus volume curve for the high pressure BCC 
structure. This is also fit with a Vinet equation of state. 
These fits are necessary to determine the phase transi¬ 
tion pressure because our current technique does not cal¬ 
culate the pressure independently. Utilizing these fits, 
we determine the enthalpy to be equal in the HCP and 
BCC phases at 807 CPa, a pressure far in excess of that 
calculated using density functional theory techniques. 

Civen the incredible sensitivity of this phase transi¬ 
tion on any errors in the calculations, it is prudent to 
investigate these errors directly. The two remaining sig¬ 
nificant sources of error are the fixed node approxima¬ 
tion for the wavefunction and the pseudopotential error 
for the Hamiltonian. In the case of solid beryllium, we 
believe that the nodal surface should be relatively sim¬ 
ple and so delay this investigation. The pseudopotential 
approximation is the more likely source of error and we 
investigate its impact by performing calculations with a 
bare Coulomb potential for the beryllium ions. In this 
sense we are following in the footsteps of Esler et ah^ 
who showed that by removing the pseudopotential ap¬ 
proximation, highly accurate results could be obtained 
for the equation of state of cubic boron nitride. 

There is a key difference between our approach and 
that taken by Esler et al. That paper used a DFT tech¬ 
nique, LAPW, to generate trial wavefunctions that na¬ 
tively satisfy the cusp condition and thereby avoid a di¬ 
vergence of the DMC local energy near each nucleus. For 
reasons of computational efficiency, we instead have cre¬ 
ated a very hard pseudopotential for beryllium that has 
all 4 electrons in the valence. This pseudopotential has 
only a single projector and a very small cutoff radius, re- 
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suiting in a potential that is very similar to the coulomb 
potential further than 0.4 bohr from the nucleus. There 
will naturally be a difference between the resulting wave- 
functions produced in the plane wave basis set used by 
quantum espresso^ and those from an all electron code, 
however we can mask this by using the freedom afforded 
by the one body Jastrow factor to reintroduce the cusp in 
the wavefunction at each nucleus. Finally, for the QMC 
calculations, we replace this pseudopotential with a -4/r 
Coulomb potential so that the resulting calculations will 
be entirely free of the pseudopotential approximation. 

This all electron methodology poses two dangers. 
Firstly, if the trial wavefunction is too poor, the timestep 
error will be incredibly difficult to converge and secondly, 
the fixed node error may be significantly increased by in¬ 
cluding the two core electrons in the valence and having 
to generate a nodal surface including them. In practice, 
neither of these concerns appear to be a problem. Using 
a timestep of 0.005 a.u. (a factor of two smaller than in 
the previous calculations) we were able to converge the 
energy to within 5 meV per atom. This is a factor of 
five larger than the pseudopotential calculations, but the 
error was consistent across all densities and phases, thus 
enabling accurate calculations of the phase transition. 
Concerns about the nodal surface are also unfounded. 
Analyzing the wavefunctions within DFT, we expect that 
the core electrons will not be greatly perturbed by the 
pressures explored in this study and therefore this error 
should be relatively constant across all of our calcula¬ 
tions. Also, the percentage of the total energy recovered 
by the DMC compared to the VMC is small. In fact it is 
on the order of 0.2% in the all electron calculations com¬ 
pared to 0.5% in the pseudopotential calculations. While 
strictly speaking this is a measure of the quality of the 
trial wavefunction globally, this is in practice highly cor¬ 
related with the quality of the nodal surface. 

The final concern in performing these all electron cal¬ 
culations is their significantly greater computational cost. 
In practice, the computational cost of the stochastic pro¬ 
cess grows only linearly with the number of electrons in 
this regime, suggesting that the dominant computation 
is the evaluation of the single particle orbitals. However, 
the variance of the all electron calculations is a factor of 
10 larger than the pseudopotential calculations, render¬ 
ing the cost of obtaining a given statistical error 200 times 
larger than for the pseudopotential calculations. Due to 
this limitation, we were unable to perform calculations 
on as large of a supercell, potentially resulting in larger 
two body finite size effects. Following Esler et ali^, we 
utilize the information on the finite size scaling from the 
pseudopotential calculations to minimize these finite size 
effects, a prescription tested by performing all electron 
calculations for the largest supercell size for one density 
in both the HOP and BCC phases. Specihcally we use 
the following formula 

^ae ~ [^ps - ^ps)^ ^ae U) 

for the HCP phase where the subscripts ae and ps refer to 


QMC energy vs volume for Be 



FIG. 3. (color online) Energy versus volume of BCC and HCP 
beryllium, including quasiharmonic phonon zero point energy. 
Insets show the pressure vs volume and enthalpy vs pressure. 
The fact that the fitting lines are on top of one another is a 
consequence of the nearly degenerate enthalpy of the phases. 


all electron and pseudopotential calculations respectively 
and the superscripts refer to the number of atoms in those 
supercells. The term in parenthesis is the estimate of the 
residual hnite size error due to performing the all electron 
calculations with smaller supercells. For the BCC phase, 
the all electron calculations utilized 54 atom supercells 
while the pseudopotential ones used 128. 

Calculating the HCP energy vs volume curve, we find a 
marked improvement in accuracy. The equilibrium vol¬ 
ume is 54.87 ± 0.04 bohr^, which is very close to the 
experimental 54.776. The bulk modulus is similarly im¬ 
proved at 115.7 ± 1.0 GPa, in excellent agreement with 
the experimental 116.8 GPa. Continuing the calculation 
to the BCC phase, we find the phase transition pressure 
is 423 GPa, nearly a factor of two smaller than the pseu¬ 
dopotential calculation and in relatively good agreement 
with DFT calculations. These results are shown in Fig[3 
The improvement in accuracy for the ambient proper¬ 
ties of the HGP phase supports the conjecture that the 
predominant error in these calculations was due to the 
use of these particular LDA derived pseudopotentials. 
The pertinent question is how to avoid this loss of accu¬ 
racy in future calculations. The idea of turning to all elec¬ 
tron calculations is appealing, but the ^5.5-6.5 gQ^ling of 
the computational cost of all electron calculations^!^ 
and the vastly increased memory required to store the 
all electron wavefunctions render this approach untenable 
for all but the lightest of elements. Glearly, better pseu¬ 
dopotentials are required. Several previous works have 
developed pseudopotentials using a different schemes. 
Predominantly there are the Hartree-Fock derived pseu¬ 
dopotentials of Trail and Needs^ Burkatzki, Filippi and 
Dolg (BFD)^ and many body correlated electron pseu¬ 
dopotentials of Trail and Needsi^ 

Here we have elected to test the BFD pseudopotential 
because it is derived from an entirely different philosophy 
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LDA PPs 

BFD PPs 

all electron 

experiment 

HCP Vo (bohr’^) 

52.27 ±0.02 

55.19 ±0.01 

54.87 ± 0.04 

54.776 

HCP Bo (GPa) 

124.21 ± 0.74 

112.99 ±0.43 

115.69 ± 1.04 

116.8 

HCP ambient C/A 

1.569 ± 0.004 

1.569 ± 0.004 

1.569 ± 0.004 

1.568 

HCP B' 

3.47 ±0.03 

3.39 ± 0.02 

3.41 ±0.06 

- 

BCC Vo (bohr’^) 

51.81 ±0.03 

54.60 ± 0.01 

54.04 ± 0.07 

- 

BCC Bo (GPa) 

BCC B' 

124.33 ± 0.21 

111.83 ±0.21 

117.61 ±0.66 

- 

3.50 ±0.02 

3.56 ±0.01 

3.60 ±0.05 

- 

Transition Pressure (GPa) 

807 

490 

423 

- 

AE (meV / atom) 

115.5 ± 1.0 

114.5 ±0.7 

99.9 ±2.7 

- 


TABLE I. Parameters for HCP and BCC beryllium calculated with LDA pseu¬ 
dopotentials, BED pseudopotentials and without pseudopotentials. Also listed are 
the calculated phase transition pressures and the change in energy per atom from 
HCP at its equilibrium volume to BCC at its equilibrium volume. 


than the DFT one, that is it is based on reproducing to¬ 
tal energies rather than the shape of the wavefunction for 
the atom. Using this potential for the current problem 
did require some additional effort as these potentials are 
presented only in semi-local form. Since we used QUAN¬ 
TUM ESPRESSO to generate the trial wavefunction and it 
uses the vastly more computationally efficient Kleinman- 
Bylander— form, we were obliged to generate projectors 
for each of the angular momentum channels. While this 
proved straightforward for Be, this is in principle a deli¬ 
cate operation as there exists a strong possibility of intro¬ 
ducing ghost states for heavier elements. The other pri¬ 
mary complication is that the cutoff radius for the beryl¬ 
lium potential (determined as the radius at which the 
potential for all angular momentum channels is within 
10“^ Ha of the bare coulomb potential) is very large, 
2.28 bohr. This means that for all but the lowest densi¬ 
ties considered, the cores will overlap. While this would 
be considered a fatal flaw for the DFT style pseudopo¬ 
tentials, calculations with BFD pseudopotentials often 
have some small overlap of the cores, while still perform¬ 
ing well. For this reason, we ignore the core overlap and 
proceed with the calculations as implemented. 

Calculating the energy versus volume relation for BCC 
and HCP using the BFD pseudopotential produces re¬ 
sults in much better agreement with the all electron cal¬ 
culations as can be seen in Table ID This is encouraging 
as it shows that it is possible to produce more accurate 
pseudopotentials for use in quantum Monte Carlo cal¬ 
culations of beryllium. However, we stop short of rec¬ 
ommending BFD pseudopotentials for all DMC calcula¬ 
tions. Three potential issues with the current methodol¬ 
ogy for applying DMC to condensed phases remain to be 
resolved. The first is that almost all modern plane wave 
DFT codes require the use of Kleinman-Bylander form, 
raising the issue of calculating accurate projectors for the 
DFT calculations needed to generate trial wavefunctions. 
The second issue is due to the potential incompatibility 
between the treatment of core electrons implied by these 
potentials and the density functional calculations that 
are generally used to generate the trial wavefunctions in 
these systems. This raises the possibility of a gross in¬ 


compatibility between the trial wavefunctions that would 
used in the DMC calculation and the pseudopotentials. 
In principle this can be resolved by optimizing all parts 
of the trial wavefunction within variational or diffusion 
Monte Carlo, however this is not currently practical for 
calculations of condensed matter. The final stumbling 
block is the large cutoff radii of the BFD pseudopoten¬ 
tials. While ignoring this overlap produced reasonable re¬ 
sults in practice for beryllium, whether this will be more 
generally true should be further investigated. 

In conclusion, we have presented highly accurate all 
electron DMC calculations of the HCP to BCC phase 
transition for beryllium, a property requiring relative er¬ 
rors to be controlled to well beyond chemical accuracy. 
These calculations found excellent agreement with the 
experimentally measured properties of the HCP phase. 
Although diamond anvil cell measurements have been un¬ 
able to detect this transition to pressures of 190 GPa, our 
calculated value is in much closer agreement with that 
obtained using density functional theory. We also found 
that the primary error in calculating the phase transi¬ 
tion pressure using the methodology of Shulenburger and 
Mattssoni is due to the pseudopotentials used. As this 
methodology is common in the QMC literature, this re¬ 
sult is of crucial importance to the field. Furthermore, we 
showed that these errors were greatly reduced when using 
BFD pseudopotentials and suggested several properties 
to be addressed in making wide use of pseudopotentials 
generated in this way for calculations of condensed mat¬ 
ter. 
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